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1 Introduction 



Baryon number is not a conserved quantity in the standard model. Rather, because of 
the anomaly, its violation is related to the electromagnetic field strength of the weak 
SU(2) group [§, 

« = Na-^e^TcF^F^ = N G £- 2 E?Bt , (1.1) 

where N G = 3 is the number of generations.^ In vacuum the efficiency of baryon 
number violation through this mechanism is totally negligible [|l| , but at a sufficiently 
high temperature this is no longer true |2|, B|. This can have very interesting cosmo- 
logical significance, since it complicates GUT baryogenesis mechanisms and opens the 
possibility of baryogenesis from electroweak physics alone. This motivates a careful 
investigation of baryon number violation in the standard model at high temperatures. 

The baryon number violation rate relevant in cosmological settings can be related by 
a fluctuation dissipation relation J|, [|, |6| to the "Minkowski topological susceptibility" 
of the electroweak theory, also called the "sphaleron rate," 



r = Jd 3 x jT eft {[E?Bf(x,t)] [EjBjM]) , (1.2) 

where () means an expectation value with respect to the equilibrium thermal density 
matrix. Here t is Minkowski time. This quantity is not simply related to the Eu- 
clidean topological susceptibility |7|], and we do not possess either perturbative tools 
or Euclidean tools to carry out its calculation. 

It has been argued by Grigoriev and Rubakov || that the value of the susceptibility 
T in the quantum theory will be the same as its value in classical Yang-Mills field 
theory. This would open a new avenue for measuring T, since classical Yang-Mills 
theory can be put on the lattice ||. There has been some progress on measuring T 
on the lattice JlO, 11], [12], |li| ; in particular two different methods have been developed 



for dealing with the right hand side of Eq. (|1.2| ) in a topological way which eliminates 
lattice artifacts in its measurement |b|, |I5f . 

At the same time our qualitative understanding of Grigoriev and Rubakov's claim 
has improved. A complication with their proposal is that 3+1 dimensional classical 
Yang-Mills theory contains ultraviolet (UV) divergences, which Bodeker, McLerran, 
and Smilga have argued may be important in setting T [[TEfl. Subsequently, Arnold, 
Son, and Yaffe have demonstrated that a particular class of diagrams, the hard thermal 



loops (HTL's), are essential to establishing T |17|1 . The amplitude of the HTL's in the 



classical theory is linearly divergent, and therefore linearly cutoff dependent. In the 



4 There is also a contribution from the hypercharge fields, but it will not be relevant here because 
the topological structure of the abelian vacuum does not permit a permanent baryon number change. 
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full quantum theory the HTL's are finite, with almost all of the contribution arising 
from excitations with momentum k ~ ttT; such high k excitations are not properly 
described by the classical theory. Arnold, Son, and Yaffe argue that because of the 
HTL's, the effective infrared (IR) theory "feels" the "cutoff" which quantum mechanics 
provides for the classical theory, and that the value of T scales inversely with the cutoff 
momentum scale. As a result, rather than the naive dimensional estimate of T ~ a 4 T 4 , 
the parametric behavior of T should be T ~ a 5 T 4 ; and in particular T is inversely 
proportional to the strength of the HTL's, which is conveniently parameterized by the 
Debye mass squared m^. On the lattice this means that T should depend on the lattice 
spacing a as T oc ag 2 Ta 4 T A , a behavior which has recently been verified numerically 



18]. Arnold, Son, and Yaffe's argument has been carefully re-analyzed by Bodeker, 
who has shown that there is an additional, logarithmic dependence on the Debye mass, 
and that, permitting an expansion in log(l/g) 3> 1, the leading behavior is actually 
r ~ a 5 log(l/^)T 4 |9|. 



If we take the limit log(l/g) ^> 1, Bodeker presents an effective theory for evaluating 
the coefficient of the a 5 log(l/g)T 4 law |l!|. The effective theory is UV safe [20] and 
the coefficient can be found accurately by lattice means []21j] . However, in practice the 
expansion in log(l/g) 1 turns out to be very poorly behaved. To get a reasonably 
accurate value for T at the physical value for the electroweak coupling, a ~ 1/30, it 
is necessary to treat the dynamics of the classical field theory with a full inclusion of 



the HTL effects. This is challenging, because the HTL effective action is nonlocal [22 



However, it is possible to rewrite the HTL action in terms of a local theory with added 
degrees of freedom, as we will discuss below. Thus, it could be possible to determine 
T by measuring the topological susceptibility of lattice regulated, classical Yang-Mills 
theory, supplemented by added degrees of freedom which correctly generate the hard 
thermal loop effects. Doing so would both test Arnold, Son, and Yaffe's claim, and 
determine the numerical coefficient of the a 5 T 4 law, and therefore tell us how efficiently 
baryon number is violated at high temperatures. 

One way of realizing this goal was presented in |23| and implemented and used to 



measure T in |24|. The purpose of this paper is to present an alternative and in some 



respects more efficient implementation of classical Yang-Mills theory plus hard thermal 
loops, and to use it to check the results of ||24j| . Our approach is based on a way of 
writing the hard thermal loops in terms of auxiliary fields which was first proposed 
in [25|. Using this formulation to incorporate the HTL action on the lattice has been 



advocated by Bodeker, McLerran, and Smilga [T(J. This paper represents a concrete 
numerical realization of that idea. 

In Section ^] we review the local formulation of classical Yang-Mills field theory 
supplemented by the HTL action due to Blaizot and Iancu and due to Nair. Their 
theory contains an infinite set of fields, so in Section |3] we perform a transformation and 
a truncation to make the number of fields in the model finite, without losing spherical 
symmetry. The resulting theory does not quite give the correct HTL equations of 
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motion; we study the difference, and how it vanishes in the limit as the truncation 
leaves in more and more fields, in Section Then we discretize space and time in 
Section ||, and review how to measure T topologically in Section ||. We study the 
numerical behavior of T as a function of the strength of the HTL's and the truncation 
point in Section [7|. 

Our conclusions are in Section [|, but we summarize them here. The HTL effective 
theory shows a dependence on the strength of HTL's which is consistent with Arnold, 
Son, and Yaffe's arguments, and grossly inconsistent with HTL independence. The 
dependence on the truncation point is surprisingly weak, so only a few new fields need 
to be added to approximate the correct HTL behavior. Thus, our algorithm proves 
quite an efficient way of incorporating HTL's. Our final results for T are consistent 
with those of Moore, Hu, and Miiller j24jl , and for the physical value of 7715 in the 
minimal standard model, = (ll/6)g 2 T 2 and a = 1/30, they give approximately 
r = (25.4 ± 2.0) a 5 T 4 . 



2 Hard thermal loops in the continuum 

In this section we discuss the origin of the hard thermal loops in terms of kinetic 
theory, and we present a local theory in which extra degrees of freedom generate the 
hard thermal loops. Nothing in this section is original; rather it is a review of Blaizot 
and Iancu's and of Nair's work |25|, |26], 27[ We include it for completeness and 
because our numerical implementation of classical Yang-Mills theory with hard thermal 
loops will be built directly from it. 

Two controlled approximations make the dynamics of IR fields in the electroweak 
theory tractable numerically, and both arise because the theory is weakly coupled. 
First, the IR degrees of freedom can to a good approximation be treated as classical 
fields. Using this fact to perform calculations of nonperturbative IR correlators was 
first proposed by Grigoriev and Rubakov [0] , and the accuracy of the approximation has 
been addressed in |L6|, [29|, j30| . The conclusion of @] is that the classical approximation 



is an excellent approximation in the infrared, but UV divergences in the classical theory 
are potentially dangerous and must be handled carefully. 

The solution to this problem is to regulate the classical theory in some way, which 
for the moment we will not specify, and then to treat the UV degrees of freedom 
separately by perturbation theory. Here the other controlled approximation enters; 
the UV degrees of freedom are described by linearized kinetic theory, up to corrections 
subleading in g. 

Since the equilibrium distribution of UV modes, iV (k), is color neutral, it does not 
directly enter in the field equations of the classical IR fields. Rather, it is necessary 
to expand the UV mode distribution function (one particle density matrix) up to first 
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order in fluctuations from equilibrium, 



N(x, k) = jVo(k) + 5N singlct (x, k) + 5N adi .(x, k) + . . . . (2.1) 

Fields in a representation higher than fundamental lead, in addition to the singlet and 
adjoint representation terms we have written, to higher representation departures from 
equilibrium; but neither these, nor the singlet deviation from equilibrium 8N singlet , 
directly interact with the IR classical fields, and at the linearized level they can be 
dropped; only No and SN^. will be relevant. Note also that N should have a spin 
index, and if there are scalar or fermionic degrees of freedom then it also has a species 
index. At leading order, corresponding to the HTL approximation, the contribution 
from each spin and species are of the same form except in the statistics for N , so we 
will not write them in what follows. 

At leading order in the coupling the IR classical fields evolve under the Yang-Mills 
field equations with a source arising from the UV modes, [25] 



{D»F V ,T = (2-2) 
= 2gC A J-0- B v lx 5N a (x,k), (2.3) 

with f M = (1, v), v = k/|k| the (ultrarelativistic) 3-velocity of the particles (note that 
f M is not a Lorentz covariant quantity), and Ca = 2 for SU(2) gauge theory. We have 
only written the contribution of gauge excitations here, there are additional terms of 
the same form for scalars and fermions where appropriate. The distribution function 
evolves via a convective covariant derivative equation which reflects the ultrarelativistic 
propagation of the UV degrees of freedom. The interactions between SN a and the IR 
classical field strength is subdominant because the coupling is weak; however, the 
electric field polarizes the equilibrium distribution, providing a source term for 5N a . 
The equation for the evolution of 5N a , at leading order in g, is 

AJ)]\I a f)]\f 

— = (v,DZy b 5N b (x, k) + gv.F^x)-^ = . (2.4) 

Note that this equation is not Lorentz covariant; it involves only the electric field, not 
the magnetic field. The reason is that the equilibrium distribution Nq has a rest frame. 
A magnetic field in that frame changes trajectories of individual particles, but it does 
not disturb the (rotationally symmetric) equilibrium distribution, whereas an electric 
field polarizes the plasma. 

One approach to making a numerical model for the IR classical fields plus UV modes 
is to simulate the distribution function N with a large number of charged particle 
degrees of freedom. In the limit that the number of particles is large and their charges 
are small, one recovers the above equations. This is the approach proposed by |23| 
and implemented in p4] . Here we will deal instead with the distribution functions. 
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This complementary approach can test the reliability of the results of |24| and may 
also prove simpler and more efficient. This is particularly true because Eqs. ( |2.3| ) and 
Q2.4D carry extra redundant information; 5N a is actually a function of k times a fixed 
function of |k|, namely, |26], 



5N a (x,k) = - g ^W a (x,v) , (2.5) 

where v = k takes on values over the unit sphere. In terms of W, the convective 
evolution of the departure from equilibrium is 

(y lt D> i ) ab W b (x, v) = v^F^(x) , (2.6) 

and the current felt by the IR classical fields is 

dtt 



J» = < J -^v,W a (x,v), (2.7) 

where dfl v means that v is integrated over the unit sphere with its natural measure. 
Here is the square of the Debye mass. These equations can be viewed as generalized 
Hamilton- Jacobi equations arising from the conserved energy^ 

H = Jd 3 x {^F^ + iifcift + \ml J ^W a (x, v)W a (x, v)) , (2.8) 

and rather nontrivial Lie-Poisson brackets [ 28 1 . Here Q v is the integration measure for 
integrating v over the sphere. 

When there are more than one species, W a represents the deviation from equilibrium 
felt by each, and is a sum of a contribution from each species of charge carrier, 

9 9 9 fN N s N f \ , s 

with N s the number of fundamental representation, complex scalars and Nf the number 
of fundamental representation, chiral fermions. In the SU(2) weak sector of the minimal 
standard model, N = 2, N s = 1, and N f = 12, so = (ll/6)g 2 T 2 . This is also a 
lower bound for all extensions of the standard model. 



Our approach will be to find a discrete implementation of Eqs. ( j2.3|) , (|2.6| ), and 
( p.7| ), and to study their evolution to determine the diffusion constant for Chern-Simons 
number. 



5 Throughout this paper Roman direction indices run over the 3 spatial directions with positive 
metric, while Greek direction indices run over all 4 spacetime indices with signature (H ). 
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3 Expansion in Spherical Harmonics 



Unfortunately, the representation of the hard thermal loops in terms of W a (x, v) does 
not provide a set of equations which are easy to implement numerically. The problem 
is that W a (x, v) is a function not only of space-time, but also over the sphere. Even if 
we discretize space onto a lattice, W still "lives" on a sphere at each lattice point, so 
it still takes an infinite amount of information to specify W completely. It is necessary 
to define W over the sphere in some way requiring only a finite number of degrees of 
freedom. Since we want to recover spherical symmetry on scales long compared to our 
lattice spacing, we should choose to do so in a spherically symmetric way. Our choice 
is to expand W in spherical harmonics, 

oo I 

W a ( X ,v)=Y: E Wim(^m(v), (3.1) 
1=0 m=-l 

where W^ m {x) is a function over space-time only. Because W a (x, v) is real valued, the 
Wi m satisfy the relations 

W l a m = (-l) m W l a ;_ m , (3.2) 

so only the real part of Wio, and the real and imaginary parts of Wi m for m > 0, 
should be viewed as independent variables. Here we use the Condon- Shortley phase 
convention and normalize Y\ m so that 

J dQYl*m Y l'm> = Sl,l>S m ,m> ■ (3-3) 

Inserting the expansion (|3.1|) into Eq. ( [2.6D , multiplying by Y{ m , and integrating over 
angles, gives the equation of motion for Wf m , 

dW a 

= -C lmjL . m .,i (A) o6 W^ m , + 5 hl v mi E? , (3.4) 
where v m i is the vector v expressed in spherical components 



dn v Y{ m (v) Vi (3.5) 

and Ci m j' m ' t i is an integral over 3 spherical harmonics: 

Q m ,i'm',i = J dVlvY^iy^iYvm'iy) . (3.6) 

We give explicit expressions for v mi and Ci m y m t^ in Appendix [A[ 
Furthermore, in terms of the spherical components the current is 

J? = ^^Wim , Jo = ^=^0 , (3.7) 
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and the conserved energy density is 



H 




-a 1 2 
Im I 



) 



(3.8) 



As written, Eqs. Q, Q, and (Q are equivalent to Eqs. O, (§, and ( P|) . 
They still contain an infinite number of degrees of freedom. However, they are in a 
form more amenable to a spherically symmetric truncation. 

The meaning of m^ ) W a (x,v) is that it represents the net charge of all excitations 
moving in the v direction at point x. What we have done is to transform to angular 
moments. rn^W^x) is the total charge of all excitations at site x. m^W^x) is 
roughly the net charge moving in the +z direction minus charge moving in the — z 
direction, and fn^Wf^ with I > 2 represent higher tensor moments in the distribution 
of excitations. For instance, a positive rn^W^ means, roughly, that there are more 
charges of type a moving either up or down the z axis than in the x, y plane. Note 
that only Wq Q and W? m interact with the IR fields, and only W% m is directly sourced 
by those fields. All the higher moments are important only in propagating the charge 
distribution through the convective derivative term in Eq. ( |3.4| ). 

The model still contains a countably infinite number of degrees of freedom, namely 
Wi m , I — 1, 2, 3, . . . . Most of these degrees of freedom are describing extremely subtle, 
high tensor fluctuations in the distribution of moving charges. It is reasonable to think 
that smearing the angular resolution of the distribution of charges by truncating the 
series of Y tm at some finite Z max will not significantly change the physics. In particular, 
for / max > 1 it will not change the way the charge current interacts with the Yang-Mills 
fields, but only the way the charges propagate; and for sufficiently large / max we expect 
the effect of angular smearing to be unimportant. Therefore, to render the set of fields 
finite, we truncate the series of W^ m at some finite Z max . The evolution equation for W^ m 
is still Eq. (|3.4j) , but with all Wy m , with /' > / max fixed to zero. Equivalently, we could 
set all Cim t i' m i t i with either / > Z max or /' > Z max to zero. The number of independent 
adjoint Wi m matrices is (7 max + l) 2 . 

As long as Ci m ,vm',i satisfies the relation 



and the terms involving v m i are either both present or both absent (they are absent if 
Zmax = 0), then the Hamiltonian, Eq. ( |3.8| ), and the phase space measure are conserved 
by the evolution equations. Hence it makes sense to speak of equal and unequal time, 
equilibrium thermal correlation functions. When / max is finite we are no longer con- 
sidering a theory which is strictly equivalent to classical Yang-Mills field theory with 
added hard thermal loops, but the behavior should approach the correct behavior in 
the limit / max — > oo and we can consider taking this limit numerically. 



lm,l'm',i — ^l' 



I'm' ,lm,i 



(3.9) 
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4 Propagator and Thermodynamics at Finite / 



Before moving on to the numerical implementation of the effective theory described 
in the last section in discrete space, we should pause to see how well or how badly 
the theory with finite / max cutoff reproduces the hard thermal loops. To do so we first 
look at whether it reproduces them correctly at the thermodynamic level; it does so 
perfectly for all Z max > 0. Then we examine the propagator of the theory, which will 
only be reproduced properly in the / max — > oo limit. 



4.1 Thermodynamics 

As discussed at the end of the last section, the theory with an / max cutoff possesses well 
defined thermodynamics described by a Hamiltonian which is quadratic and diagonal 
in the W/ m 's. The only complication is that the phase space is constrained due to 
Gauss' law: 

(D.Ey = ^Lw o a , (4.1) 



and the partition function reads 

Z = J VA{DE{DW lm 5 ((£> • E) a - m 2 D W^/V^) exp(-if/T) , (4.2) 
H = \j* x (b»B« + E1E1 + ^ E |W7J 2 ) • (4.3) 

Every W\ m except W^oo is Gaussian and they can all be integrated out immediately. It 
is also convenient to introduce a Lagrange multiplier for Gauss' law, 

5 ((£> • E) a - mlW^I VA?j = J VA exp {iA a [{D ■ Ef - m 2 D W^ /V^} /T} . 

(4.4) 

Doing so makes E and W^o Gaussian as well, and they can now be integrated out, 
yielding 

Z = JvAiVA ex.p(-H'/T), (4.5) 

H' = yd 3 x[BtBt + (D t A r(D t A r + m 2 D A a A-], (4.6) 

where the wave function term for Aq arises from integrating out the E field and the 
Debye mass squared term arises from integrating out W^o- 

The sole thermodynamic consequence of the W fields is the introduction of a Debye 
mass, and its magnitude is given exactly by the coefficient in the W field equations of 
motion. This corresponds exactly with what the complete hard thermal loop thermo- 
dynamic contribution should be. Furthermore, the Debye mass is introduced even for 
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Imax = 0, the absolute minimum value. We do not recommend using Z max = 0, how- 
ever, because in this case the W fields have no dynamics and every Wq Q is a conserved 
quantity. Therefore, the system is not ergodic and a Hamiltonian trajectory will not 
densely sample the microcanonical ensemble. However, to the best of our knowledge 
the only conserved quantities (besides the Gauss constraints) in the non-abelian theory 
with Z max > 1 are energy and momentum^, and we expect ergodicity in this case. The 
conclusion is that the technique reproduces the thermodynamics of the full HTL theory 
exactly, for all Z max > 1. 

4.2 Propagator 

Now we turn to the study of the propagator in the Z max cut-off theory. We work only to 
linear order, or equivalently, we will study the propagator only in the abelian theory. 
In this case we can study one k mode in isolation. 

Since the spherical harmonic expansion does not break rotational invariance (even 
when we restrict I < / max ), it is sufficient to study the propagation of modes for which 
k is in the 3-direction. The Fourier transformed equations of motion are 

2 

u 2 A m -k 2 A m = ~f-W lm (4.7) 
wWim - kCi mi i> m > t3 Wi> m > = uSi tl A m , (4.8) 

where we have defined A m= ±i = ^ An /Q(^Ai + iA 2 ) and A m=0 = ^JAtt/3A 3 . The A±i 
are the transverse components of the gauge field and A m=0 is longitudinal. 

Since C[ m ^ m i ^ oc S m>rn i, c.f. Eq. ( |A.4j ), the equations of motion do not mix different 
m-sectors (this is the advantage of choosing k || e 3 ). We also note that W/ max j max 
and Wi nmx -i nmx do not evolve at all. In general, the components with m ^ ±1 do 
not couple to the transverse gauge fields. We will not be concerned here with the 
propagator in the longitudinal sector, or with any sector which does not couple to any 
gauge fields, so the only 'interesting' modes are those with m — ±1. It should be noted 
that this decoupling occurs only in the abelian theory. (It also allows a more efficient 
representation for hard thermal loops than the one we use here, see [ftH|.) 

In the following we choose m = 1, which is the sector which couples to the transverse 
gauge fields. The matrix Cw = Cn,/'i,3 is a symmetric and traceless matrix of size /^ ax 
with non-zero (positive) elements only if I' — I ± 1. (Note that, because |m| < /, / is 
restricted here to the interval 1 < I < Z max , hence the dimensionality of Cw-) As a 
result, in the eigenvalue problem 

C X a = * a X a , (4-9) 

the eigenvalues X a are real and non-degenerate, and they come in positive and negative 
pairs: if A is an eigenvalue, so is —A. If Z max is odd, the matrix has one zero eigenvalue, 

6 On a discrete lattice the total momentum is not conserved, due to the Umklapp-effect. 
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otherwise the eigenvalues are non-zero. The eigenvectors x a are rea l an d orthogonal, 
and we will normalize them to be orthonormal. 

Writing the matrix Cw in terms of the eigenvectors and eigenvalues, 



we can solve for Wn in Eq. 



C«' = ExrA°X?, (4.10) 



Inserting this in Eq. (|4.7|) , we obtain the inverse transverse propagator 



A-i = _ u 2 + e + Hh±J2—^(x1) 2 . (4.12) 

(max 3 _ fc^O! ^ ' V ' 

Q = l 



Let us now compare the propagator (|4.12 ) to the theory without the Z-cutoff. Re- 
member that in this case the propagator has a cut in the interval — k < uj < k [p2 |, and, 
in the limit uj <C A; <C m^, it describes overdamped behavior with damping coefficient 
k 3 jm 2 ^. Thus, the damping rate is ~ g A T when k ~ g 2 T, which is the relevant 
momentum scale for non-perturbative physics. 

What does the propagator look like at different values of Z max ? If = 0, the gauge 
fields are decoupled from the W fields, except through Gauss' law (see Eqs. ( |2.3|) and 
(|3.7|)); transverse physics is the same as in the absence of the W fields. At / max = 1 the 
"matrix" Cw = is a scalar, and the propagator describes a massive vector particle: 
Af 1 = -uj 2 + k 2 + m 2 D /3. 

The first interesting case is Z max = 2. The propagator is still easy to solve analytically, 
and (using Eq. ( |A.4| )) the inverse propagator becomes 

The propagator has two zeroes given by uj 2 = k 2 /5, and 4 poles at 



, 3k 2 ml 1 {6k 2 ml\ 4k 4 

uj 2 = 1 -±- A 1 . (4.14) 

5 6 2 \ ^ 5 3 ) 5 



In the limit k 2 <C m\ the poles are 

^ = !^l + h 2 + 0(k'), u 2 = ^ + 0(k e ). (4.15) 
3 5 5 mf) 

The first 2 poles correspond to the plasmon, and give it the right dispersion relation 
up to corrections of order k^/m 2 ^. The second 2 poles are at uj ~ g 3 T for k ~ g 2 T and 
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Poles of propagator l m 



-0.5 




Figure 1: Left: The inverse propagator Eq. (|4.12| ) with Z max = 10, plotted against 
w/mo with fixed k = OAmi). Right: The positive frequency poles of the propagator 
at / max = 10. In these figures, one can clearly see the development of the cut in the 
interval —k < uj < k, and the two plasmon poles at uo 2 pa m 2 D /3 + 6/c 2 /5. 



m D ~ gT. Thus, instead of the correct overdamped behavior, the Z max = 2 propagator 
( |4.13| ) describes oscillatory behavior with u ~ g 3 T. At first sight, this may look like 
a fatal flaw in the /-mode cutoff method. However, as we will argue below, in practice 
this is not a serious drawback. 

For odd values of / max the matrix C\y has one eigenvalue equal to zero. As with 
/max = 1, the self-energy contribution to Eq. ( |4.12| ) has a constant 'mass term', 

^(Xi 0) ) 2 + ^ E — J TtM?- (4-16) 

What this means is that there is a linear combination of W and A fields, namely 
Wn = Wxf\ A = {m 2 :) /?)k 2 )Wxi \ which is strictly static. Thus, part of the "power" 
in the A fields is lost to the dynamics of the system. There are also propagating modes, 
both at the plasmon frequency and for to < k. For Z max = 3, the poles are at 

u 2 = H± + Q k 2 + 0{k 4 h u 2 = 8_ k 2 + 0{k 4 ) _ (41?) 



We can identify the same plasmon pole as with Z max = 2, Eq. Q4.15| ), but the other 



pole behaves as |a>| ~ k instead of |co>| ~ k 2 . For relevant values of k, the poles of the 
^max = 2 propagator are at much smaller \u\ than for Z max = 3. 
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This pattern is seen to be true also for larger / max . While we have been unable to 
find a general analytic expression for the poles of the propagator, it is easy enough to 
solve the eigenvalue problem ( [4.9D and find the poles of the propagators numerically. 
In Fig. [I] we show the inverse propagator and the location of the poles when Z max = 10. 
In general, we can state the following about the poles of the propagator: 

(i) For Z max even, there are Z max poles and / max zeroes of the propagator in the interval 
—k < uj < k. For odd values of l max , the number of poles and zeroes is / max — 1. In 
either case, as Z max — > oo, the poles and zeroes merge into a cut in the propagator. 

(ii) There is a pair of plasmon poles at uo 2 Pd + (6/5)k 2 + 0{k A /m 2 ^). 

(iii) When Z max is even and k <C mc / v^max, the lowest pair of poles behaves as 

uo^± (4.18) 

m DV'max 

whereas the other poles in the region \u\ < k depend linearly on k. For Z max odd, 
all of the poles in this region are linear. As we make Z max larger, the power lost 
to the static mode becomes smaller roughly as / max . 

The absence of cuts means that the gauge field propagation is non-dissipative. We 
should expect this behavior in the abelian theory because the equations are linear. 
However it need not concern us, because at large Z max the behavior differs from the 
^max = oo limit only over very long time scales, and the nonlinearities in the non- 
abelian case should become important on shorter time scales if / max is sufficiently large. 

The spectral power density 

p(u,k)/u; = (2/w)ImA(a; + ie,k) (4.19) 

for fixed k = 0.4mD is plotted in Fig. 0, both for the full propagator without the Z-cutoff 
and for several (even) values of Z max . In the finite Z max case the spectral density gets 
contributions only from the poles of the propagator ( 4.12|) : 



2tt 

k )/u = 6 ( u - u vo\e{k)) x Res A imax (w pole (/c), k) . (4.20) 

poles 

The spectral power is strongly concentrated around ui = with a peak width 5uj w 
4k 3 J '(Trmjy) . The spectral power of the Z max propagator closely follows the / max = oo 
curve; however, in order to have enough power in the central peak region, Z max should 
be large enough so that there are poles well within the bulk of the peak, which is the 
relevant region for the propagator to describe the correct damping. 

We can use this property to derive an approximate 'rule-of-thumb', which tells how 
large Z max should be for a given value of k (which, for the relevant physics, should be set 
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k/m n = 0.4 





o.o 

-0.40 -0.20 0.00 0.20 0.40 

G)/m D 

Figure 2: The spectral density p/ui at k = 0.4mo for propagators without the /-cutoff 
and with various values of Z max . The spectral density for finite / max is a sum of form 
J2 a C a 5(ui — uj a ). The plot symbols are plotted at coordinates (u a , 2C a / {u a +i —uj <x -i))i 
which makes it possible to compare the different Z max -values. 

to g 2 T): we simply require that Z max is large enough so that lowest positive frequency 
pole satisfies 

4k 3 

oo P oic(k) < 5- . (4.21) 

Numerically, this corresponds to 

C x n > 0.62m 2 D /k 2 - 0.8, > IMml/k 2 - 1.1 , (4.22) 

with good accuracy. Strikingly, one has to use 3 times larger values for Z max in the 
odd sector than in the even one. This is due to the lack of the u ~ k 2 -pole in the 
odd sector, as emphasized above. While for modest values of ^ 2 T/m D = 0.3. . .0.5, 
^max = 6 or 4 should be sufficient, for very weak coupling or large m-o the required 
^max-value becomes impractical for numerical work, since the numerical effort will rise 
as (I 

max 

Naturally, one has to remember that Eq. (|4.22 ) is based on an ad hoc requirement 



that the lowest positive frequency pole should be within the peak of the spectral den- 
sity, and different criteria would lead to very different requirements (but the overall 
pattern in Eq. ( [4.22|) should remain) . What value of / max one really needs in non-abelian 
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simulations may differ from Eq. (|4.22|) by a large factor; indeed, the numerical results 
for non-abelian theory in Sect. [7] seem to imply that Eq. ( }4.22| ) is overly strict. 

We also note that the poles of the W field coincide with the A field poles, so we have 
gotten them for free. For odd / max , there is one pole not accounted for yet: the mode 
corresponding to the Cw zero eigenvalue, which does not propagate at all. Naturally, 
the spectral power of the W propagator is very different from the A propagator. 

So, what do these results tell us about the sphaleron rate in the non-abelian theory? 
We conclude the following: 

A. When Z max is odd, a component of the gauge field is static and not fluctuating, 
and therefore does not contribute to real time processes. Since the static component 
is largest in the infrared, we expect this to reduce T relative to large / max limit. This 
behavior is worse for small / max and should go away at large Z max as the static component 
contains less and less of the total gauge field amplitude. 

B. In the even / max sector, the location and density of the poles is relevant for 
the correct damping in hot plasma: the larger Z max is, the more the poles are able 
to reproduce the concentration of spectral power at small u, and so the stronger the 
damping and the smaller the sphaleron rate. Thus, when Z max increases, the sphaleron 
rate should approach the physical one from above. The approach should be much faster 
than in the odd Z max sector, see Eq. (|4.22|) . 

This behavior is indeed close to what we observe in Sect. [7|. 

Obviously, the results in this section imply that for fixed Z max one cannot have the 
correct leading order (in g) behavior of the gauge field propagator in the strict small 
g limit. We expect this to be true also for the non-abelian gauge propagator, and 
hence for the sphaleron rate. Nevertheless, for realistic values of g and mo we expect 
a modest / max to be sufficient. 



5 Lattice equations of motion 

In this section we discuss the discretization of the continuum equations of motion 
Eqs. (|2.3|), (|3.4j ) and (|3.7|). Naturally, not all of the properties of the continuum evo- 
lution can be satisfied on a discrete lattice, but the update rule of the lattice system 
should fulfill at least the following criteria: 

(i) Gauge invariance, lattice translational and rotational symmetry and C, P, and T 
symmetries are preserved, 

(ii) Gauss' law is identically satisfied, 

(iii) The total energy is conserved. 

Naturally, we also require that the small lattice spacing and smooth field limit gives 
the correct continuum behavior. 

The discretization of the system is very similar to the pure Yang-Mills theory, de- 



veloped by Kogut and Susskind [33]. The lattice is a 3-dimensional torus of size 
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L 3 = N 3 a 3 , with lattice spacing a. As is customary in real-time simulations, we 
use A = gaugeQ and discretize the gauge fields in terms of spatial parallel trans- 
porters Ui(x) = exp(igaAi) e SU(2), and electric fields Ei(x) which belong to the Lie 
algebra of SU(2). Ui(x) and Ei(x) live on the links connecting points x and x + i (here 
we use the shorthand x + i for x + ae^). The Wi m fields are located on lattice sites. 
Thus, for each lattice site the total number of field variables is 3 SU(2)-matrices and 
3 + (/max + l) 2 adjoint matrices. 

On the lattice we want to use dimensionless field variables. We absorb the lattice 
spacing and g in lattice fields as follows: 

gaA — > A , ga 2 E — > E , gaW^W. (5.1) 

For compactness, we also use dimensionless lattice coordinates, integer, 
reintroducing a when necessary. We shall consider the evolution of the lattice fields 
both in continuous and discrete time. In discrete time, one update step consists of 
evolving the fields from time t to t + S t , where 5 t <C 1 in order to keep the evolution 
stable and integration errors small. 

5.1 Gauge field update 

We shall use the standard single plaquette definition for the magnetic field strength: 

l r l 



- J <t\r -{■■;■{■■•■ - (3 h £ [1 - |Tr Uu] ■ (5.2) 



Here U a is the ordered product of the link variables around a plaquette, 

Ua,ij(x) = Ui(x)U s (x + i)Uj(x + j)Uj(x) , (5.3) 

At tree level /?l = 4/(^ 2 Ta). However, this receives radiative corrections; these will be 
discussed below. 

Let us now consider the lattice gauge field equations of motion both in continuous 
and discrete time. The (continuous) time derivative of the link matrix U is given in 
terms of the electric field as0 

d t Ui{x, t) = iEi{x, t)Ui{x, t) . (5.4) 



7 We emphasize that this choice is just a convenient way to fix the gauge ambiguity in the field 
update laws, and that any alternative choice would give the same value for gauge invariant correlators. 



Ei(x) appears on the left in Eq. (5.4) because we choose to record Ei(x) so that it transforms 
under gauge fields as an adjoint object at the basepoint x rather than the endpoint x + i of the link 
from x to x + i. Alternately we could work in terms of Ei{x) = Uj(x)Ei(x)Ui(x), in which case 
the expression would involve UE rather than EU; similar changes would appear in other expressions 
involving E. There is no physical difference between the two choices. 
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The 'gauge force' term in the evolution equation of the electric field is fixed by the 
magnetic Hamiltonian ( |5.2|) , by varying Eq. ( p.2| ) with respect to Ai. When we add the 
current term due to the W\ m fields, we obtain the evolution equation for Ef 



d t E?(x,t) = -z|Tr r a U t (x } t) J2 Sjj&t) 
1 \j\¥* 



+ \mx,t)+V?jt(x + i ,t)]. (5.5) 



Here Sij is the gauge link 'staple' which connects the points x and x + % around the 
plaquette: 

S tj = UjWUiix + fiUfo + i) (5.6) 

The summation index j in Eq. ( |5.9|) goes over both positive and negative directions; a 
negative value means that the link is traversed in the opposite direction as in Eq. (|5.6j ): 

U ,(•'•) /•;(•'" ./)• 

The current terms in Eq. (|5.5|) are given by jf = (mDa) 2 /(47r)f^ i H / 1 a m . Since Ei(x) 
is located between the points x and x + i, the current ji(x) is averaged between the 
beginning and the end of the link. The current at x + i has to be parallel transported 
to point x, and we use the shorthand expression 

Vf <$>\x + i,t) = [Ui(x, t)$(x + i, t)Uj(x, t)] a (5.7) 

for the adjoint field parallel transport from point x + i to point x.0 

In discrete time, the adjoint field Ei transports the link matrix Ui(t) to Ui(t + 5t). 
In order to keep the evolution symmetric in time, it is natural to place in the half- 
timestep value t + ^5t- Integrating Eq. ( |5.4|) , we obtain the discrete time evolution 
equation for Uf. 

Ui{x, t + 5 t )= exp [iEi(x, t + \8 t ) S t ] Ui(x, t) . (5.8) 

Alternatively, one can think of exp(iEi5 t ) as being the timelike plaquette in the (t,i) 
plane, which updates U as shown because we have chosen A Q = gauge. (In another 
gauge there would be an extra A dependent term in the U field update, and in the 
updates of the E and W fields as well; it is the convenience of leaving these out which 
encourages the choice of temporal gauge.) 

The discrete time electric field update can be obtained now from Eq. ( |5.5| ) by sub- 
stituting 

dtE?(x,t) - yiEt(x,t + \S t ) - E?(x,t- §<5 t )] . (5.9) 



The lhs of Eq. ( p.5|) remains as is even at discrete time. As formulated, the discrete 
time update steps (|5.8|) and ( |5.9|) are symmetric under time reversal and they give an 
algorithm accurate to order 0(5f). 



9 If we worked in terms of E, the other current would require parallel transportation to the end 
point of the link. 
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As mentioned above, the relation /?l = 4/(g 2 Ta) receives corrections because UV 
modes behave differently on the lattice than in the continuum. This has been calculated 
in Appendix B of |24j (see also f34|), with the result 



0i 



1 37A U 2m 2 D a 2 m> 4 \ £(m D a) 



g 2 aT \3 6vr J \3 3 18 J 4vr 

i 2 D t 



1 m 2 D a 2 \ S(mDa) 



\3 18 J 4tt 
Here £ = 0.152859 . . ., and S(moa) and ((wiDa) are integral functions: 



(5.10) 



S(m ) = r d3k 1 — = r d3fe A 1 — , (5.H) 

1 ; J-n (2vr) 3 k 2 + m 2 ' ^ ^ i-w (27r) 3 (P + m 2)2 1 > 

where k 2 = X^4sin 2 fc»/2. To 5% accuracy, this can be expressed as /?l — 4/((7 2 aT) + 
0.61 for values of mDfl used in this work. However, we shall use the full expression in 
our analysis. In the sequel we will write /3l for the variable appearing in Eq. ( |5.10| ), 
and write (3 for 4/g 2 aT. 

Further subtleties related to this thermodynamic correction arise when we convert 
T to continuum limits; we will address this in Appendix |C|. 

5.2 Wi m update and doublers 



The Wi m equation of motion Eq. (|3.4|) has only first order derivatives in time and 
space. In order to preserve the exact P and T symmetries on the lattice, the first 
order derivative terms should be replaced by symmetric finite differences. Thus, the 
continuous time lattice equation of motion for Wi m is 

d t Wi m (x,t) = ~ Ci m>Vmlii [ViW Vml {x + i,t) - V-iW Vml (x - i,t)] 

+ ^ 1 v mi [E t (x,t)+V^E i (x-i,t)]. (5.12) 

The electric field contribution is symmetrized from each of the links which connect to 
point x. 

As was done with the spatial derivative, we substitute the time derivative dtW with 
a symmetric finite difference [W(t + S t ) — W(t — S t )]/ (25 t ), and the value of W at time 
t + S t will depend on values at times t and t — 8 t . Explicitely, the update rule becomes 
a 'leapfrog' 



Wi m (x, t + 5 t ) = Wi m (x, t-5 t ) + 5 t ^2,5i^v mi E mCti 

- C^^^ViW^ix + i^-V-iW^ix-i^Y (5.13) 
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Here E ave is the average electric field influencing the propagation of W\ m from t — S t 
to t + 5 t . Since this is over two timesteps, there are 4 timelike 'plaquettes' to each 
direction i: 

E aYe>i (x,t) = ^[Ei(x,t - \5 t ) +V-iEi(x -i,t - \S t ) 

+ Ei(x, t + \S t ) + V-iEi(x -i,t+ \5 t )} . (5.14) 

Note that, due to Eq. fl5.8|) , the parallel transport V-iEi(x,t + \5 t ) can be made with 
U matrices either at time t or time t + St with the same result. In practice, one does 
the Ei transport once for each timestep, and stores the result for the next timestep. 
To summarize, the discrete time update step (t — > t + 5 t ) goes as follows: 

1. start with U{t), E(t - S t /2), W{t) and W{t - S t ), 

2. evaluate E(t + 5 t /2) with Eqs. Q and Q, 

3. calculate W{t + 8 t ) with Eq. (^JJ) (and forget W{t - 5 t ) and E(t - 5 t /2)), and 
finally 

4. calculate U(t + S t ) with Eq. 

A generic feature of a first order differential operator on a discrete lattice is the 
decoupling of 'odd' and 'even' coordinate sectors: Wi m (x,t + 5 t ) depends only on Wi m 
at points (x, t—S t ) and (x±i, 5 t ); in particular it does not depend on Wi m (x, t), which is 
its immediate predecessor. More precisely, if we label the coordinates with an integer 
valued parity label p = J2i x i + t/5 t , the Wi m fields at odd and even values of p do 
not interact, except through their coupling to the gauge fields. This causes a species 
doubling problem, in analogy to the one familiar from lattice QCD (the Dirac equation 
is of first order). The properties of the doublers in a linearized theory are discussed in 
detail in Appendix p. 

There are 15 extra low-energy doubler modes, living around the corners of the 4- 
momentum space hypercube hi = (0,7c/a), uo = (0, 7r/(5 t a)), with at least one of k iy uj 
non-zero. The continuous time equation of motion ( |5.12|) has only 7 spatial doublers; 
the rest are introduced by the time discretization (and can be avoided, see subsection 
|5.4j) . However, in contrast to lattice QCD, in our case the doublers are benign: first, 
they couple only very weakly to the gauge fields, decoupling completely at the corners 
of the Brillouin zone (see Appendix ||). Second, they couple only to gauge fields at 
very high wave numbers k ~ 1/a and/or frequencies uj ~ l/(a5 t ). Thus, the doublers 
do not influence at all the physically interesting small k and uj gauge field dynamics, 
and their effect on modes close to the lattice cutoff remains small. 

Because the time step is small (5t 1), the timelike doubler modes uj ~ ir/aSt are 
especially weakly coupled to low-frequency modes. Indeed, in simulations we used S t = 
0.05 and observed no appreciable energy transfer between the timelike doublers and 
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low-frequency modes. However, since the timelike doublers are low-energy excitations 
of W fields which are not present in continuous time, they can cause problems in 
thermalization of the system and, as it turns out, in counting the active degrees of 
freedom. This will be discussed below in section |576\ 



Before leaving the update we should comment on energy conservation. In continuous 
time, we can write down the lattice version of the Hamiltonian (51!): 



H{t) = - |Tr [/□(*)] + WE?(x,t) + £ \W lm (x,t)\ 2 . (5.15) 



This Hamiltonian is exactly conserved by the equations of motion (571), (|5.5| ) and 



( |5.12| ). However, in discrete time there is no equivalent conserved expression. A good 
approximation to the energy can be obtained by symmetrizing the contribution of the 
electric fields in Eq. ( |5.15| ) with respect to t: 

E?(x, t) - [Ef{x, t - \5 t ) + Ef(x, t + \5 t )\/2 . (5.16) 

The energy obtained this way fluctuates with an amplitude oc 5 2 , but the mean value is 
stable. The conservation of mean energy is guaranteed by the time reversal symmetry 
of the discrete time equations of motion: if, at some point in the evolution of the fields, 
we invert the sign of E (E — > —E) and conjugate and reverse sign for the hard particle 
charges (Wi m — > —W^ = — (— l) l Wi m ), the system will exactly retrace its evolution 
backwards. If the energy had a tendency to increase, inverting the time would cause it 
to decrease. Since the configurations (U, E, W) and (U, —E, — W*) are just as likely to 
appear in a thermal distribution, the system cannot exhibit any systematic tendency 
for the average energy to change. The stability of the system is a necessary property 
for long Hamiltonian evolutions. 



5.3 Gauss' constraint 

The Gauss' law is given by the 0-component of the equations of the motion 



D t F*°=j 



H^oo- 



(5.17) 



On the discrete spatial lattice and discrete time, care has to be taken to make the 
appropriate symmetrizations to the fields F i0 = Ei and W appearing in Eq. (|5.17|) . 
Since E^ is living on half timestep time values t + |5 t , we symmetrize H^o from times 
t and t + 8 t : 



Ei(x, t + \8 t ) - V-iE^x -i,t + \8 t ) 

i 

(mDfl) 2 1 



+ 



4tt 2 



[W 00 (x,t)+W 00 (x,t + 5 t )] = 



(5.18) 
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This condition (or rather, the constancy of the violation of this condition) is satisfied 
exactly by the evolution equations (|5.8|), ( |5.9| ) and ( |5.13| ). To see this consider the 



change of Eq. (|5.18|) under one time step. It gets contributions from each dE/dt and 
from dWoo/dt. (There are no contributions from the time derivative dll/dt of the U 
appearing in the parallel transporter V-i because dll/dt commutes with E and cancels 
between the U and U> in Eq. ( |5.7| ).) In the absence of W fields, the time derivative of 
Eq. ( |5.18|) is zero, as shown by Ambj0rn and Krasnitz |10|. The addition of W fields 



adds new terms to the Wqo field and E field updates. First there is a contribution to 
Ei(x) and V-%Ei(x — i) from W\ m {x). According to Eq. ( |5.5|) these are equal; but Ei(x) 
and V—iEi(x — i) appear in Eq. (|5.18 ) with opposite sign, so there is no contribution 



here. There is also no contribution to dW Q0 (x)/dt due to Wi m (x). Second, W\ m at 
each neighboring site contributes both to dE/dt on the link between the neighboring 
site and x, and to dWw/dt, through Eqs. ( |5.9|) and ( |5.13| ) respectively; but the two 
contributions to the time derivative of Eq. ( |5.18| ) cancel, because Coo.im.j = v* mi . Hence 
the update preserves Gauss' law if it is satisfied by the initial conditions. Enforcement 
of Gauss' law is therefore a problem for the thermalization algorithm, not the evolution. 

5.4 A way to eliminate temporal doublers 

There is an alternative way to write the update rules which eliminates all the high 
frequency doubler modes, which we now discuss. First, note that the reason there 
are doublers is that the update as specified in the previous subsections requires and 
maintains twice as much information about the W fields as is necessary. As discussed 
in the summary at the end of subsection [B~2l , the update needs the value of Wi m at two 
time slices. However, only W and not its time derivative appear in the Hamiltonian, 
so a complete specification of the fields should only require Wi m to be specified once 
at each site. The excess information describes the state of the doublers. Eliminating 
the doublers will require eliminating half of this information. This is possible since, 
as noted earlier, the update rule for W does not mix the W fields on odd and even 
sublattices. Therefore, it is possible to define W only at every other spacetime point; 
we can define it only at the even sites, that is, points for which p = [t/5 t + J^i^i] is 
even. Eq. ( |5.14j ) remains unchanged, but Eq. (|5.9| ) has the modification 

ji(x,t) , t/5 t + J2iXieven 

[j?( x ,t)+V? b j>(x + i,t)]^{ , (5.19) 

V?fl(x + i,t) , t/St + EiXi odd 



that is, we use whichever j is defined. Similarly, in Gauss' law, Eq. ( |5.18| ) involves 
either Woq(x, t) or Wqo(x, t+St), whichever is defined. The time derivative of the Gauss 
constraint remains conserved, for the same reasons as before. 

Updating the fields in this way removes 8 of the 15 doublers and cuts the number 
of computations, and hence the CPU time, almost in half. It may slightly increase 
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timestep errors because of the even-odd alternation of the current in the E field update 
rule; but this can be compensated for by reducing 5 t , which is not problematic because 
of the reduction in the number of computations per time step. We have compared 
the update with and without this modification and find that the results for physical 
measurables agree within statistical errors. 



5.5 Lattice thermodynamics 



In continuous time the equations of motion ( |5.4|) , (5.5), and (|5.12|) describe a Hamil- 
tonian evolution which conserves energy and phase space volume. We can study the 



thermodynamics of the system by using the Hamiltonian (|5.15| ) to write down the 
canonical partition function 

Z = J IjldUi^dEiix)] [ J] dW lm (x)} l[S(G(x))e- H/T , (5.20) 



xlm 



where G(x) is Gauss' law, Eq. (|5.18|) (in continuous time). Introducing a Lagrange 
multiplier field Aq in exact analogy with what we did in continuous space in subsection 
TO, we can integrate out the E and W fields to obtain the lattice partition function 



Z 



[ \li dU i( x )] \li d M 

J x,i L x 

/3 L £[l-|Trf/ D (t)] + 



x 



-Ha 



(5.21) 



(5.22) 



The gradient term for the A field is the simplest lattice implementation of the contin- 
uum (DiA ) 2 , and appears as the A mass term without any corrections, just as 
in the continuum case. The form of the partition function above is equivalent to the 
path integral of the full quantum theory in the high-temperature dimensional reduction 
approximation on the lattice |35], j36[] . This guarantees that this theory reproduces the 
(equal time) thermodynamics of the Yang-Mills fields. 

This property can be used to fix the bare lattice value of the mass term mo- In 
general, classical field theories suffer from UV divergences; however, when we consider 
the static thermodynamics of the theory in Eq. fl5.22| ), only a finite number of UV 
divergent diagrams appears. These divergences can be absorbed in counterterms, and 
in particular for the theory in Eq. ( |5.22| ), we have p5 



ni 



D,bare 



m 



D,phys 



E 5 2 T 



na 



3.17591... 



(5.23) 



Here ?«D,phys i s fixed according to the actual particle content of the theory, see Eq. (|2] 
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5.6 Thermalization 



The real time simulation has to be started from a configuration which has been chosen 
from a thermal distribution so that the the Gauss' constraint is satisfied. As emphasized 
above, to start the update we need the fields U(t), E(t — 5 t /2), W{t) and W(t — 5 t ). 

We will use the same general philosophy as in fi~H| . Some of the degrees of freedom, 
namely Ei and Wi m , are Gaussian, while others, namely Ui, are not. We can draw 
the Gaussian fields from the thermal ensemble and then use the evolution equations to 
"mix" this thermalization with those degrees of freedom which are not Gaussian. The 
thermalization proceeds by evolving the Hamiltonian equations of motion of the system, 
but periodically "refreshing" the Gaussian degrees of freedom, that is, discarding the 
values of Gaussian degrees of freedom and drawing them from the thermal ensemble. 

At first sight, this plan appears to be complicated due to the Gauss' constraint. In 



the case without the W fields this problem was solved in |LT[ , by first drawing E from 
the Gaussian distribution ignoring the constraint and then projecting to the constraint 
surface. It is trivial to extend that technique to the current situation. However it 
is actually possible to do something even easier. Only the component Wqq of Wi m 
enters the Gauss' constraint. Thus, according to Eq. ( |5.15|) , we can set the higher Im- 
components freely to the correct thermal distribution, that is, draw each of W^, I > 1, 



from a Gaussian distribution of width ^ 8tt / (m^a) 2 . The thermalization then proceeds 
as follows: 

(1) Set U(x,t) = 1, E(x,t) = Woo{x,t) = 0. 



(2) Choose Wi m {x,t), I > 1, from the Gaussian distribution of width J 8n / (m^a) 2 . 

(3) Evolve the equations of motion for a short period, transferring energy from W\ m 
to the other fields, while preserving Gauss' law. 

(4) Repeat from (2) until the fields are thermalized. 

However, in discrete time we do not have an exact Hamiltonian, and there is an 
inherent ambiguity oc 5 2 in the definition of energy. It is not immediately evident how 
the fields should be thermalized. At a more practical level, the randomization of W\ m as 
above is complicated by the fact that we need Wi m fields at times t and t— St to start the 
leapfrog update. This is closely associated with the timelike doublers of the W fields. 
However, as was discussed in section [5.2|, the timelike doublers couple extremely weakly 



to the low-frequency mode sector, and there is practically no energy transfer between 
the two sectors. This was also seen in simulations: the energy contained in the doubler 
modes remained at the level where it was set by the initial thermalization during the 
whole trajectory.^ Moreover, the gauge fields care only about the low frequency modes 

10 More precisely, the energy transfer remains negligible for timestep St = 0.05 used in the simula- 
tions. Using a dangerously large time step of order St ~ 0.2-0.3, energy transfer becomes significant. 
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(see Appendix ||). Thus, in principle, we are at liberty to do whatever we choose about 
the timelike doubler modes; we can either thermalize them or try not to excite them in 
thermalization. The gauge fields will not see the difference — however, in the former 
case W fields will contain roughly twice as much energy as in the latter. Note also that 



the whole problem would go away if we used the update discussed in subsection 5.4 



In all of our 'production' runs we chose not to excite the timelike doubler modes. This 
makes the lattice modes resemble as closely as possible the continuous time fields. Note 
that the Hamiltonian (|5.15| ) counts the degrees of freedom and energy equipartition 
correctly only if there are no timelike doublers, and the 0(5%) ambiguity in energy is 
valid only in this case (in the presence of doublers, the ambiguity is of order 100%). 

The thermalization without the doublers can be accomplished using the steps (1)— (4) 
as above, but replacing the step (2), for example, by one of the following two methods: 

(a) Set Wi m (t) to Gaussian random variables in step (2) above, and perform the first 
update step in (3) using a forward asymmetric time difference for these Im modes: that 
is, instead of approximating the time derivative with [W(t + 5 t ) — W(t — 5 t )]/(25 t ) in 
Eq. ( |5.13|) , we use [W(t + 5 t ) — W(t)]/5 t . This is a natural way to start a leapfrog, 
and it gives a smooth interpolation for the fields. This method gives slightly incorrect 
mean energy, but the error is 0(5%). 

(b) Set Wi m (x,t — 5 t ) = Wi m (x,t), where I > 2, to Gaussian random variables in step 
(2). Now also the first step can be performed with the leapfrog. This method excites 
the doublers more, but the amplitude of their excitation is only 0(5%). Note that now 
only modes / > 2 can be randomized, since in one timestep both E and Woo interact 
with with I = modes. 

In our production runs we used the method (b). We also made test runs with the 
doubler modes fully exited. This is simple to accomplish: proceed as in items (l)-(4) 
above, randomizing only Wi m (t), I > 2, and perform the evolution with the leapfrog 
update ( |5.13| ). Since there are now twice as many active Wi m modes, the width of the 
Gaussian distribution has to be multiplied by v^2, in order for the U and E fields to have 
the same total energy as before. As mentioned above, in the gauge field observables 
the doublers have no observable effect. 

Let us note that a Langevin-type thermalization, as used in [57 for pure Yang-Mills 
theory, would be straightforward to implement by coupling the noise to Wi m fields. 
Indeed, coupling the noise only to the highest /-modes might be of interest even during 
a simulation, since this could mimic the effect of the higher I modes. 



6 Measuring the Chern-Simons number diffusion 

The baryon number violation rate is related to the diffusion of the Chern-Simons num- 
ber, defined as the charge associated with the right-hand side of the anomaly equation 
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(0): 

Ncs = SrJ Sx ^ k { F * Al ~ \f^ A ) AC ) = ^J d " xJ B ■ 

Since SU(2) has a non-trivial third homotopy group 7T3(SU(2)) = Z the Chern-Simons 
number Nqs is a topological index for vacuum configurations: we can perform any 
gauge transformation to a trivial configuration A = without any cost in energy, and 
the resulting configuration is as good a vacuum configuration as the initial one. iVcs 
is equal to the winding number of this gauge transformation. Since it is now integer 
valued, it classifies the vacuum configurations into disconnected classes, which cannot 
be continuously gauge transformed to each other. Thus, a vacuum-to-vacuum process 
which increases Nqs smoothly by one unit must go through a non- vacuum exited state, 
the sphaleron. Due to the axial coupling to fermionic current, this process lifts one 
left-handed solution of the Dirac operator from negative to positive energy, and pushes 
one right handed state from positive to negative energy. Since the SU(2) sector of the 
standard model is a chiral theory which does not couple to the right-handed fermions, 
this process will create one fermion for each fermionic generation (Nq). 

At high temperatures the Chern-Simons number diffuses readily, and integer values 
are not particularly preferred. In any given volume the Chern-Simons number performs 
a random walk in time, and the diffusion constant, T, can be measured from 

r = Mm lim - «!> . (6 . 2) 

Here the angle brackets (•) refer to an average over the thermal ensemble. The change 
in the Chern-Simons number can be evaluated from 

N CS {t) - iVcs(to) = #i /"* dt' I d 3 xE?B? . (6.3) 

57T Z Jto J 

In principle, this measurement is readily convertible to lattice language: lattice ver- 
sions of the fields E and B feature prominently in the equations of motion (|5.4j ),( pT5|) . 
However, this "naive" definition of Nqs on the lattice, often used in the early work on 
Chern-Simons number diffusion in lattice SU(2) gauge theory || [10], [IT], suffers 



from spurious noise and diffusion which obscures the physical iVcs diffusion. Moreover, 
due to its UV nature, the amplitude of the noise diverges as 1/a in fixed physical 
volume, which is disastrous in the continuum limit. The reason for this noise is well 
understood: the integral over lattice E ■ B on the right-hand side of Eq. (|6.3|) does 
not form a total time derivative, and hence it depends on the path along which one 
connects the initial and final configurations in Eq. (|6.3| ). In other words, it does not 
give us a topological measurement. 

In general, topology of lattice fields is ambiguous, since the variables are always 
continuously connected to trivial ones. However, at fine enough lattice spacings (still 
easy to achieve in our simulations) almost every one of the plaquettes is very close 
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to unity in a thermal ensemble; large plaquette values are exponentially suppressed. 
Perturbatively this means that the gauge fields are small, and for this subset of lattice 
fields topology can be unambiguously defined |j38| . Physically, this means that the 



spatial size of the topology changing configurations, sphalerons, is large in lattice units. 
This will be true because the energy of sphaleron-like configurations increases linearly 
with inverse size. The Boltzmann suppression factor for small (lattice scale) sphalerons 
is enormous and they, in practice, never appear in simulations. The interplay between 
entropy and the Boltzmann factor sets the dominant sphaleron size to be ~ g 2 T. For 
topology to be unambiguously defined, our lattice spacing must be considerably smaller 
than this. 

Two successful methods for measuring topology in the current "real time" context 
have been recently developed. The first method uses an auxiliary "slave field" to 
track the winding number of the gauge WM. It is a development of a method originally 



proposed by Woit [p9| , which is based on counting the winding numbers of singularities 
in the Coulomb gauge. In this work we use the second method, "calibrated cooling." 



This method is based on the cooling method by Ambj0rn and Krasnitz [13|, and fully 
developed by Moore in [I5|. The rest of this section will summarize this method. 

The calibrated cooling method relies directly on the fact that the sphalerons are 
large and extend over several lattice units. Thus, we can get rid of most of the ultra- 
violet noise in the thermal configuration by applying a small amount of cooling to the 
configuration: the resulting configurations are very smooth on lattice scales, but they 
still have the same topological content as the original configuration. After cooling the 
integral ( |6.3|) can be performed with small errors. The accumulation of residual errors 
is prevented by periodically cooling all the way to a vacuum configuration: we know 
that the true vacuum-to-vacuum SN CS = integer, and any deviation is due to accu- 
mulated integration error, which can thus be "calibrated" away. This is schematically 
described in Fig. |3|. 

The cooling path is defined by the gradient flow of the standard single plaque- 
tte Kogut-Susskind gauge action, given in Eq. ( |5.2| ). The evolution along this path 



is parametrized by fictitious cooling 
equation of motion is now JT3|] 

dUi(x) . „dA? 



"time" r, dimensionally (length) . The cooling 



dr 



la 



dr 



■Ui 



x 



10~ 



a 1 



Tr 



ia a U,-, 



X 



£4 



X 



Ui{x) 



(6.4) 



Here the staple S is defined as in Eq. ( |5.5|) . On the lattice the equation above is evolved 
in discrete r, and we use here optimized step lengths by alternating 5r/a 2 = 5/48 and 
10/48. Too large a time step causes the UV modes to become unstable. 

The evolution of (6.4) all the way to a vacuum configuration is a computationally 
demanding task, and it can easily dominate the cpu time. The integration can be 
dramatically accelerated by blocking the lattice: after a bit of cooling the fields are 
very smooth at the lattice scale, and essentially no information is lost if we reduce the 
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Hamiltonian time 



configurations generated by the Hamiltonian trajectory 




integrate Ncs along cooled trajectory 



integrate Ncs 

from cooled trajectory 

to vacuum 



~ vacuum configurations 



Figure 3: How the Nqs evolution is measured (after ||15|| ). Top horizontal line shows 
the configurations (solid circles) generated by the lattice equations of motion. Every 
few timesteps, the configurations are cooled a fixed cooling length, giving a parallel 
cooled trajectory (open circles). Now the fields are smooth enough so that E • B can 
be reliably integrated, giving SN C s(t) along the cooled trajectory. In longer intervals, 
the Nqs measurement is "grounded" by cooling all the way to a vacuum configuration. 
If SNqs along paths like V\ — > A — > B — > V% is always close to an integer, we know 
that the integration errors are small. The residual deviation from an integer value is 
subtracted from SNqs(A — > B), cancelling the accumulation of errors. 



number of lattice points in each direction by a factor of 2. The blocked [/-matrices are 
formed from the product of the matrices on the two links between the blocked points. 
Since the lattice spacing a is now increased by a factor of 2, computational cost in 
cooling is reduced by a factor 2 5 - a factor of 4 coming from the increase in the St 
step. In our calculations we block the configuration twice in the course of cooling to 
vacuum. For detailed information about this method we refer to |I5|. 



In all of our simulations we cooled a configuration from the Hamiltonian trajectory 
at intervals St = 0.5a (once in 10 timesteps with St = 0.05a). We cooled to depth 
r = a 2 45/48 for the cooled trajectory (see Fig. ^), using unblocked configurations. SN CS 
was integrated along this trajectory using improved 0(a 2 ) accurate definitions for E • B 



TTJ. The cooling to vacuum was performed with an interval St = 12.5. Our parameter 
choices were overly conservative: the vacuum-to- vacuum integration error was typically 
of order 0.02-0.04. Thus, it is possible to use much more aggressive optimization than 
we use here without losing the topological nature of this measurement (see ref. 
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7 Simulations and results 



Our aim here is to answer the following questions: 

(1) what is the dependence of the Chern-Simons number diffusion rate T on the finite 
Z max cutoff, and is there an Z max which is 'large enough' for practical purposes or 
is an / max — > oo extrapolation necessary? 

(2) is T, in physical units, independent of the lattice spacing? 

(3) how does T depend on the physical quantity mr>/g 2 T? 

Let us first discuss the relation between the physical Debye mass and the bare 
mass parameter m^a on the lattice. As explained in Sect. |5.5| , the bare mass receives 
renormalization counterterms and diverges in the UV limit as 1/a. However, according 



to the scaling arguments of Arnold, Son and Yaffe fl7 |, the sphaleron rate should not 
actually depend on the Debye mass, which characterizes static screening properties of 
the hot plasma, but on the damping rate of the transverse gauge field propagation. 
As explained in Sect, |L this is related to the Debye mass in the continuum. However, 
due to the lattice dispersion relation, the hard gauge field modes do not propagate at 
the speed of light, and their effect on the damping is reduced. Averaging over all of 
the directions of the lattice momenta, Arnold ^0] has calculated that the effect of the 
hard gauge field modes on the lattice is a factor of (0.68 ± 0.2) times smaller than the 
continuum relation between the damping coefficient and would imply. The error 
quoted is systematic, and it takes into account the rotational non-invariance of the 
lattice propagators. Thus, we shall use the following relation between the bare lattice 
tbd and the continuum one: 



ZJn^DMtt = m D, P h y s - 0.68^-^ . (7.1) 



Here is a radiative correction of form 1+0 (a), see Eq. ( |(J.14] ). We use the improved 
relation Eq. ( |5.10|) to relate the lattice spacing a to the physical scale g 2 T. There 
are additional radiative corrections associated with renormalization of the lattice time 
scale, which we discuss in Appendix |C[ In order to avoid the uncertainties associated 
with the UV counterterm, we use mostly fairly large physical values of m-Q so that the 
UV term remains subdominant. The results are actually quite robust against variations 
in the numerical coefficient 0.68, even with the smallest mo we use. 



/ max dependence: In order to study how the sphaleron rate depends on the value 
of /max, we performed a series of runs with 24 3 lattices using fixed /?l = 8.7 and 
mj=) = 1.5/a 2 = 7.9g 4 T 2 , and varied Z max from to 10, as shown in Table |l[ The 
results are also shown in Fig. |J When Z max is even, the results are remarkably stable: 
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run parameters 


Z max time/ a 


r/(« 4 T 4 ) 


3t = 8 7 Via 3 = 24 3 
mf, = 1.59/a 2 = 8.20c/ 4 T 2 


8000 

1 20000 

2 20000 

3 30000 

4 20000 

5 20000 

6 30000 
10 20000 


1.49(15) 

0.0606(70) 

0.531(34) 

345(23) 

0.520(22) 

0.445(30) 

0.534(27) 

0.518(34) 


ft = 12.7, V/a 3 = 32 3 
mf, = 1.98/a 2 = 20.7# 4 T 2 


2 37500 
4 37500 
6 45000 


0.249(23) 
0.198(18) 
0.199(18) 


(3 L = 12.7, V/a 3 = 32 3 
mf, = 0.29/a 2 = 4.84# 4 T 2 


2 37500 
4 37500 


0.839(36) 
0.687(50) 



Table 1: How the sphaleron rate F depends on I 



indeed, the data from / max = 2 to 10 are mutually compatible within the statistical 
errors. However, for odd Z max the rate remains substantially smaller, approaching the 
even sector value from below when / max increases. 

The special case Z max = has a rate which is ~ 3 times larger than the / max = 2,4,... 
rate. This is actually close to the rate measured from standard SU(2) gauge theory 
without any W fields at the same /?l UlSR ; there the rate was 1.68 ± .03. 

This behavior is qualitatively in accord with the theoretical analysis in the abelian 
theory in Sect. f|. The odd Z max sector gives a substantially reduced rate because much 
of the infrared power is in non-propagating modes, and is therefore not available to 
participate in Chern-Simons number diffusion. However, for the even Z max sector we 
do not see the gradual decrease in the rate as predicted by the analysis in Sect. |], the 
rate just snaps to the correct level immediately when the damping is turned on by 
going from Z max = to Z max = 2. According to the requirement for minimum / max given 
in Eq. fl4.22| ), we should use / max > 0.62171^/ g 4 T 2 — 0.8 ~ 7 (for even l max ). The naive 
limits given in Eq. ( |4.22j ) are obviously too strict for the non-abelian theory. 

At larger m^/g^T 2 the difference between Z max = 2 and higher values should be 
more visible. Indeed, in simulations at m^/g^T 2 = 20.1, using Z max = 2, 4 and 6, we do 
observe a significant decrease in the rate as Z max increases from 2 to 4; this is shown in 
Table [1]. Here we use a smaller lattice spacing, /3l = 12.7, and correspondingly larger 
volume in lattice units. In this case the required Z max , according to Eq. ( |4.22| ), would be 
~ 12. We also see an effect in the rate at m^/g^T 2 = 4.75, /?l = 12.7, using Z max = 2 
and 4. 
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8.2 gf p L = 8.7 V=24 3 



1.5 



i / even 

max 

) / odd 

max. 



1.0 



0.5 



2X 



0.0 



10 



Figure 4: The dependence of T on Z max on a lattice of size 24 3 ,/? L = 8.7. 



/3l 


^max 


V 


(m D a) 2 




time / a 


T/a 4 T 4 


k' 


8.0 


4 


24 3 


0.375 


2.63 


10000 


1.23(11) 


40.5(3.6) 


8.7 


4 


24 3 


0.766 


4.68 


37500 


0.808(33) 


47.5(1.9) 


8.7 


2,4,6,10 


24 3 


1.59 


8.20 


90000 


0.526(15) 


54.2(1.6) 


8.7 


4 


24 3 


3.51 


16.4 


25000 


0.230(18) 


47.5(3.7) 


12.7 


4 


32 3 


0.291 


4.84 


37500 


0.687(31) 


46.8(1.9) 


12.7 


6 


32 3 


0.707 


8.74 


20000 


0.417(48) 


45.8(5.2) 


12.7 


4,6 


32 3 


1.97 


20.7 


82500 


0.199(13) 


51.7(3.5) 



Table 2: The Chern-Simons diffusion rate V and the parameter of the Arnold-Son- 
Yaffe scaling law k', Eq. ( |7.5|) . If the results at different values of l max are statistically 
compatible, we have taken an average over them, as indicated. 



Physical sphaleron rate: The results for the sphaleron rate are shown in Table 
and plotted in Figs. | and || The "old argument" p, |J, based on dimensional analysis, 
says that the rate should scale as F = /ta 4 T 4 with k a constant. This behavior is 
clearly excluded, as already seen in [Tj|. Rather, the rate falls linearly in g 4: T 2 /m^ ) , 
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Chern-Simons diffusion 



1.50 



A(3 = 8 
•Po = 8-7 
■ P =12.7 



1.00 



8 
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0.00 




0.00 



0.10 



0.20 
(g T/m D ) 



0.30 



0.40 



Figure 5: The Chern-Simons number diffusion rate F in physical units. Dashed line: 
fit to linear + second order term, Eq. (|7.2|); continuous line: fit to linear + a log-term, 



Eq. ([fj) 



confirming the ASY scaling picture. Indeed, we can make a fit of form 

,4tt2 / ^4t-t2 \ 2 



r 



a 4 T 4 



ci — — + c 2 



mf 



ml 



(7.2) 



to the data, with the result c\ = 4.5 ± 0.2, c 2 = —3.2 ±1.1, with x 2 — 12 for 5 degrees 
of freedom. Within our statistical errors, we did not observe any systematic lattice 
spacing dependence, and we use the results obtained with all the lattices in Table |2] in 



the fit. If we include the known logarithmic contribution pi 

Tlo; 



o 4 T 2 

(0.425 ±0.027)^-^ log 



ni 



D 



ml 



4 T 2 



« 4 T 



(7.3) 



we can perform a one-parameter fit 

r s 4 t 2 



a A T 4 



ml 



0.425 log 



' m D 
4 T 2 



+ d 



(7.4) 
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Figure 6: The Chern-Simons number diffusion rate T, expressed as k' = 
(T/a 5 T 4 )(m^ ) /(g 2 T 2 )). Dashed line: fit to linear + second order term, Eq. ( |7.2| ); con- 
tinuous line: fit to linear + a log-term, Eq. (7\4). For comparison, we include the 
results obtained with the 'particles' method and also without any hard thermal 
loop degrees of freedom |fL8| . These points are not included in the fits. 



with d = 3.09 ±0.08, with y 2 = 15 for 6 degrees of freedom. The logarithmic contribu- 
tion actually makes the fit a bit worse; however, a subleading term 0(g 4 T 2 /m^) would 
not change it, since its coefficient would be compatible with zero. The errors quoted 
above are purely statistical; we shall discuss systematic errors below. 

The rate becomes approximately constant, and the difference between the two fits 
becomes more visible, if we plot the rate in terms of the coefficient of the ASY scaling 
law k,', defined through 



m 



(7.5) 



D 



The values of k' are given in Table |2| and shown in Fig. |6|. We also include here the data 
calculated by Moore, Hu and Miiller with the particle degrees of freedom inducing the 
hard thermal loop effects [24~|, and the results obtained by Moore and Rummukainen 
using only SU(2) gauge fields without any additional hard thermal loop degrees of 
freedom fll8"f. In the latter case the damping arises solely through the UV gauge field 
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modes on the lattice, and m 2 ^ can be obtained from Eq. (|7.1| ) by setting latt = 0. 

The consistency of the results obtained with different methods in Fig. [5] is remarkable. 
This gives strong credibility to the view that the damping of the sphaleron rate seen 
in the simulations is really caused by physical hard thermal loop effects. Perhaps 
surprisingly, the pure Yang-Mills results are perfectly in line (within the statistical 
errors) with the results obtained with the hard thermal loop effective theories, even 
though in the former case the spectrum of the hard modes is strongly distorted by 
the lattice dispersion relation f40fl . Also, the consistent decrease of k' with increasing 
(g 2 T/mE>) 2 in Fig. | strongly suggests that this subleading effect is not due to lattice 
effects, since the damping is due to very different mechanisms in theories with or 
without additional hard thermal loop degrees of freedom. 

As can be seen from the x 2_va hies reported above, the quality of the fits shown in 
Figs. H and |6] is poor. This is primarily due to the 'pull' of the /?l = 8.7, m 2 D = 8.2 g 4 T 2 
-point, which has very small statistical errors. Exclusion of this point would make the 
fits acceptable; however, we do not have any a priori reason for rejecting this point, so 
we keep it in the analysis. We take the badness of the fits into account by expanding 



the errors in the quantities reported below by a factor yx 2 /V, where v is the number 
of the degrees of freedom in the fit. 

In the limit m 2 D — > oo the coefficient of the ASY scaling law becomes k' = 55.9 ± 
3.5 using the polynomial function in Eq. ( |7.2| ). However, more relevant for physical 
applications is the Standard Model (MSM) value at m 2 ^ = (ll/6)g 2 T 2 and a w = 1/30. 
This corresponds to point (g 2 T /m£>) 2 = 0.23 in Fig. ||, and thus there is no need to 
extrapolate in m 2 D . 

Finally, as discussed in the beginning of this section, the numerical coefficient 0.68 in 



Eq. (7J.) has an estimated (quite conservative) systematic error bar ±0.2. This error 
has little effect on n' if we extrapolate to m 2 D — > oo, but at the physical MSM value 
it actually gives the leading contribution to the total error. When we take this into 
account, we obtain the physical value 

k'(MSM) = 46.6 ± 2.0 stat ± 3 syst , (7.6) 

and the MSM Chern-Simons diffusion constant becomes (with combined statistical and 
systematic errors) 

T = 25.4±2.0a 5 T 4 . (7.7) 
This value is in perfect agreement with the results obtained both with the particle hard 



thermal loop degrees of freedom |24j] and with the classical Yang-Mills theory Jig] . 

It is actually likely that the systematic error of the coefficient (0.68 ±0.2) in Eq. ( [7.1|) 
is overestimated: the mutual consistency of the results obtained with and without 
the extra hard thermal loop degrees of freedom becomes noticeably worse when this 
coefficient is more than ~ ±0.1 away from the central value. 
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8 Conclusions 



Classical Yang-Mills theory plus hard thermal loops is the IR effective theory for the 
SU(2) sector of the Standard Model above the electroweak phase transition, and it 
should be used to determine the "sphaleron rate" T, which sets the efficiency of baryon 
number violation. 

We have developed a numerical implementation of the method of auxiliary fields, 
originally developed by Blaizot and Iancu and by Nair J25|, ^7], ^8[. The auxiliary 
fields are expanded in spherical harmonics and the series is truncated at a finite / max ; 
then the theory is put on a lattice. The resulting numerical model is an efficient and 
systematically improvable representation of the desired effective theory. 

Within errors we observe no lattice spacing dependence, and the convergence to the 
large Z max limit is surprisingly rapid. This means that the lattice numerical model is 
both accurate and efficient. Using it, we verify the Arnold- Son- Yaffe scaling behavior 
for r 0, T = K'(g 2 T 2 /m 2 D )a 5 T\ 
(ll/6)g 2 T 2 and a w = 1/30, the rate is 



If we use the Standard Model values of 



r = (25.4±2.0)a 5 T 4 ~ (1.05 ± 0.08) x l(r 6 T 4 



(8.1) 



The final result is in good agreement with the results previously obtained by Moore, 
Hu, and Miiller |23]. It is also in agreement with the results obtained in pure lattice 
Yang-Mills theory [pl| using the matching technique developed by Arnold j|0| to relate 
r in pure classical lattice Yang-Mills theory to its value in the quantum theory. 

The problem of determining the sphaleron rate in Yang-Mills theory is settled, at 
least at the 20% level. 
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A Spherical coefficients v m i and Ci m y m i i 

In this appendix we give explicit expressions for the coefficients Ci m j> m > ti , Eq. (|3.6|), 
and v m i, Eq. ( |3.5|) . We use the conventional normalization for the spherical harmonic 
functions: 

/ dQvY^Yiw = 5u>5 m , m > . (A.l) 
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The coefficients v m i can be given simply as 



dn v Y* m (v) Vi 

I An [An _ [An 

5 m ,i + 5 m -i) + i]l -jrr8i,2{8 m ,i + S m -i) + y — 5 ii3 5 mi o . (A.2) 



The matrix elements 



Y,v Mi f ^^(v)y 1M (v)y iW (v) 

AT 



are conveniently expressed in terms of the spherical components 

Clm,l'm',l — ^lm,l'm' + ^im,i'm') 

Cl m ,l'm r ,2 = ~^[plrn,l'm' + ^lm,l'm') • 

Finally, we can write 

^hn,i'm' = Aty^m^Si-ijiSm-irf — A(l,—m)8i + i t i'8 m -i tm i 
Cim,vm' = A(l,—m)8i-i } it8 m +i,m'—A(lim)8i+ij'8 m+ i jm i 
Cim,i'm',3 — B{V ,m)6i-. lt i>8 mim i + B(l,m)8i +li i<8 mirn i , 



where the coefficients are 

A(7,m) 

B(Z,m) 



1/2 



(Z + m + l)(Z + m + 2) 
2(2/ + l)(2/ + 3) 

(/-m+l)(/ + m+l) nl/2 
(2/ + l)(2/ + 3) 



(A.3) 



(A.4) 



(A.5) 
(A.6) 



B Linearized lattice propagator 

In this appendix we shall study the properties of the linearized gauge field propagator 
with hard thermal loops on the lattice, as was done in section |] in continuum. As 
emphasized in section ||, the second-order 'leapfrog' update for W^ m decouples the even 
and odd parity sites from each other. Here parity p = Y,i x i + t/8 t . This decoupling 
creates extra low-energy poles, doublers , in the gauge field propagator. Here we show 
that these doublers are not relevant for the gauge field dynamics. 
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Following Sect. [| we linearize Eqs. (|5.8|) , (|5.9| ) and (|5.13|) , and study one k mode in 



isolation. In general, this need not be parallel to any of the major lattice axes. In 
Sect. |5| the spherical harmonics were written in "lattice basis" coordinates, that is, the 
Y[ m components were mapped to lattice coordinate axis directions in the customary 
way. Naturally, there is no fundamental reason (only great convenience) to do this, 
and here we choose to parametrize the spherical functions as in Sect. |], so the Y\ m 
"a^-direction" is parallel to k. Then the transverse fields should oscillate only along 
the plane defined by m = ±1 components. 

Let us make a Fourier transformation of the lattice equations of motion; after some 
work we obtain the equations 

2 

U 2 A m -k 2 A m = ^J= D mm> W lm > (B.l) 
W^lrn " ^iQra/m'^i'ra' = /^D mM A M . (B.2) 

Here n = k/k, and the lattice momentum functions are 

~ 2 kid „ 2 uj8 t a 

hi = - sin , us = - — sin , 

a 2 o t a 2 

ki = - sin kid, ^ = 7 — sin u5 t a, (B.3) 
a o t a 

and the matrix D mm >, m = 0, ±1, is defined as 

ka I 3 

D mm , (k) = V 7™ cos -i- 7^ 7™ = W —Rij(n)v mj (B.4) 

t 2 V 4vr 

Rij is a rotation matrix which rotates n parallel to the lattice 23-axis, and v m i is 
defined in ( |A.2| ). The matrix 7 can be understood as a transformation between the 
lattice coordinates and the n-based spherical coordinates. Here 7^7 = 77^ = 1. The 
cos(fcja/2) factor in D mm i arises from the spatial symmetrization in Eqs. ( |5.9| ) and 
( p,13| ), and without this we would have D — 1. Also, due to the timelike symmetrization 
Gj instead of uj appears on the rhs of Eq. (|B.2|) . 



Due to the matrix D mm i the equations of motion do not diagonalize to independent 
m-components, as in the continuum. However, if k is parallel to any of the lattice axes 
we have 

D mm ' = <5m,m' ( <Wl + ^m,0 COS — 1 . (B.5) 

In this case D = 1 for transverse modes. Now we can solve for the transverse gauge 
field in Eqs. ( B.l|) , ( B.2|) as in Sect. ^ and we obtain the lattice version of the inverse 
propagator in Eq. (|4.12|) : 

,2 



+ ^ 2 + ^E— ^r(ff) 2 - (B.6) 

<J a=l U) — k\ a 
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Figure 7: Left: the positive frequency poles of the lattice gauge propagator within the 
Brillouin zone, as functions of k, for mo = 0.5 /a. For each of the continuum poles 
near origin there are doubler poles at the corners of the Brillouin zone. The thick line 
is the plasmon, and the dashed lines are the timelike doublers. For clarity, this figure 
is plotted with unrealistically large 5 t = 0.75. Right: The spectral power (residue) of 
the poles when 5 t = 0.1. The lines are as in the left panel. None of the doublers carry 
a significant fraction of the gauge field propagation. 

The pole structure of the propagator is not immediately evident from Eq. ( |B.6| ). Never- 
theless, the propagator has new doubler poles for each physical low energy (uj ~ 0, k ~ 
0) pole, as shown in Fig. [7| These are located in the corners of the momentum square 
< ka < 7r, < uj5ta < n. 

In Fig. [7] we also show the spectral power 

p(u, k) = Im A(w + ie, k) (B.7) 

of the poles (residue of the propagator). At momenta close to the lattice cutoff n/a 
almost all of the power is carried by the plasmon pole, which does not have doublers. 
Thus, at high momenta the gauge field essentially decouples from W fields, except for a 
mass term which equals mo/v^ for k along a lattice axis (but not everywhere, it is zero 
at ak = (tt,tt,tt)). This also occurs in the continuum. The poles other than plasmon 
are significant only around the physical k, u ~ corner. Interestingly, even here the 
plasmon has a power which is a factor of ~ 5 larger than the other poles. However, it 
is these poles which are significant for the non-perturbative small-frequency physics. 

One might worry about the small-/c temporal doublers, which correspond to modes 
which flip sign at each consecutive timestep: after all, these can have arbitarily long 
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spatial wavelength. These are shown in Fig. [7| with dashed lines. However, it turns 
out that these poles have even much less power than the spatial doublers. Thus, the 
existence of the temporal poles should not affect the gauge field behavior at all. Indeed, 
even in non-abelian theory simulations, where the fields are fully interacting, we saw 
no significant energy transfer between the temporal pole sector and the 'normal' small 
frequency sector. 

The propagator becomes much more complicated to study when we do not require 
that k is along any of the lattice axes. However, the pole structure of the propagator 
is qualitatively similar to any direction, and it will pick the full complement of poles 
at each corner of the Brillouin zone. 

In order to make the left panel of Fig. [7] readable, it is plotted using 5 t = 0.75. 
This brings the frequency spread of the poles to the same order of magnitude than the 
separation between the temporal doublers and the other poles. For a more realistic 
S t < 0.1 the poles would lie almost along ua5 t /ii = 0, 1 lines. 



C O(a) matching for T 

In this appendix we compute 0(a) radiative corrections which arise in the infrared 
dynamics of the lattice theory due to the compact nature of the gauge action and 
the manner in which the original equations were discretized. The goal is to find what 
modifications must be made to T and (where here is really being used to 
represent the magnitude of the damping rate for gauge field modes with u <C k ~ g 2 T, 
which determines the relevant dynamics |17], ; when we write we mean the value 



which gives the same damping rate using the continuum relation between and the 
damping rate). 

C.l Corrections to t and 

To begin, we discuss the relation between the lattice and continuum values for E, DjFji 
(meaning the first term on the right hand side in Eq. ( |5.5| )), and time t. Where possible 
we will suppress spatial and group indices, in particular we refer to DjFji as DF. We 
will write E^ etc. for the lattice fields scaled to continuum units directly using Eq. 
( |5.1|) , but always using a as given in Eq. ( |5.10|) . The scaling between the continuum 
A field and the lattice one, defined as U = exp(igaA a T a ), is gauge dependent, and we 



will always use the continuum normalization. The calculation here relies both on |34 
and on Appendix A of |ST| very heavily. 



Define the following renormalization constants: 

Z„ = /3//3 L ~ 1-0.61//?, (C.l) 

^ + 6 ^)1^ 1 + - 314/ ' 3 - (&2 > 
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N S 

Zw = 1 -— — - 1 - 0.2527//3, (C.3) 

where Z 9 is computed in |34| and presented above in Eq. ( |5.10|) , Ze is first computed 
in the appendix of [ETJ (where it has the unfortunate notation of (1 + corr) 2 ), and Zy/ 
is new to this paper and discussed more in the next subsection. 

To begin observe that, just from its appearance in the Hamiltonian next to we 
have 

(El) = Z g (E 2 c ) =► E L = Z^E a . (C.4) 
Also, from Appendix A, we have 



^ = 4 /2 ^ = Z^Z^Eo =► t L = Z E ^Z^H c . (C.5) 

We apply this correction when we extract the continuum value of V from data which 
appear as a time series in t L , so T quoted in this paper is always scaled by Vt c the 
continuum volume and time. Next, to find the renormalization of DF, we can use pure 
gauge theory relations 

d -§ L = ~DF L and = -DF C , (C.6) 

at l ate 

to find that 

DF L = Z E /2 Z g DF Cl (C.7) 

which we will need below. 

To compute the radiative corrections in the W field contribution to the gauge field 
damping rate we need to consider the equations of motion of the full system. If there 
were only infrared fields, then the first errors from our discretization (sampling neigh- 
bors to determine a derivative , . . .) would enter at 0(a 2 ), while here we will only 
be interested in 0(a) effects. Nevertheless, because of the different behavior of UV 
modes on the lattice than in the continuum, three new corrections arise: one in the 
E field source in the W equation of motion, Eq. (|5.12|) ; one in the W field source in 
the Yang-Mills-Maxwell- Ampere equation, Eq. (|5.5|) ; and one in the W field convec- 
tive covariant derivative in Eq. ( |5.12| ). The former two occur because, whereas in the 
continuum these equations relate fields at the same point (see Eqs. ( ft .4]) and (|3.7Q ), on 
the lattice they involve averages over nearby points, see Fig. || The latter correction 
occurs because the derivative term necessarily involves the gauge field connection. In 
each case the gauge field connection enters, and the interaction receives tadpole con- 
tributions which are absent in the continuum. As a result, the effective IR equations 
of motion look (in a simplified notation, dropping all subscripts including I, m indices 
and all Clebsch-Gordan coefficients, which are the same as in the text of the paper) 

= -DFL-ml^WL, (C.8) 

atj, 
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Figure 8: 
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Neighbor averaging involved in the updates of E due to W and W due to 



dW r 



k 2 E l - k 3 v ■ D C W L 



dt L 

KiK 2 = Z W , K 3 = Z^r 



(C.9) 
(CIO) 



Here is the lattice parameter converted to physical units by scaling with factors 
of a. We derive the size of the corrections k 1j2 3 in the next subsection. 
Re-arranging Eq. (|C.9|) to 



K *dT L +v - D ° 



W = Kn 1 K 2 E L 



(C.ll) 



formally inverting it, and substituting the solution for W into Eq. ( |C.8| ), gives 



dE L 
~dT T , 



-DFt 



-id ^ 

" 3 dF L +v - Dc 



KlK 2 K? 1 E L . 



(C.12) 



It has been argued in |4T|, ^] that in the overdamped case it is permissible to drop both 
the dE/dt term and the time derivative appearing in the inverse operator. Technically 
doing so commits an error of order 0{g i T 2 /m^). Note however that errors of precisely 
this size already arise from subleading corrections to the hard classical lattice mode 
contribution in Eq. ( |7. 1|) . Therefore, in Figure || there is an unknown systematic error 
in the slope of the fit line, which we will not be able to eliminate. However, we can 
still ask to make all 0(a) corrections which would affect the intercept. To do so we are 
permitted to drop the time derivatives mentioned above, giving 



/2 



[v ■ D 



.i dA 



c 



c 



dt 



DF, 



c 



c 



(C.13) 
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which gives us the 0(a) renormalization appropriate for the W field Debye mass term, 
namely, the value to use as the strength of the gauge field damping term is 

Z'lml = ZtfzgPz-W ml . (C.14) 
C.2 Evaluating ^12,3 

We see from Fig. [8] that the product arises from the difference between W a (x) 
and {l/4)[PfW b (x + ia) + 2W a (x) + Vl b W b (x - ia)]. We compute this difference for 
a very slowly varying W field in Coulomb gauge. (In this gauge the effects of the A 
field on the dynamics do not differ between the lattice and continuum, see pi| .) 
The parallel transport V° b W b [x + id) is 

T c (U t (x)W + uKx)) c = (l + igaT a A?(x)-^T a T b A?(x)A b (x) + ..) 

xT c W c x j^same, i -iV (C.15) 

and on expanding (and writing (l/A)(W(x+ia)+2W(x)+W(x — ia)) as W) eventually 
gives 

Vf W b (x + ia) + 2W a (x) + V ab W b (x - ia) 



W *T a + ^W b (4 c (x) - A\{x - ia)) f abc T a 
9 2 a 2 



+ y —W a (A\{x)At{x) + A b (x - ia)A c i (x - ia)) [T c , [T a ,T b ]] + 0(a 3 )(.C16) 

In momentum space the first term here is —f abc (ga 2 /4) J l liW b (l — k)A-(— I). It cannot 
lead to a contribution proportional to W(k) because {f abc W a (k)W b (l - k)A c (-l)) = 0. 
Therefore the first term does not rescale the interaction, so it does not contribute to 
k,ik,2- However, the last term does lead to renormalization of W a T a , of magnitude 

(faT f d 3 {ak) 1 / kj\ 2 

fadbfbdeT — J [l ~ ^ j COS (afc/2) , (C.17) 

where k = (2/a) sin(a/c/2) and the integral runs over the Brillioun zone, aki e [— tt, 7r]. 
Using identities from [^|, the value of the integral is (1/2)S/(4tt). Therefore the final 
rescaling we find is 

KxKi = Zw = \-^. (C.18) 

The calculation of k 3 proceeds similarly. Here we need to compute (Vf'W (x + 
ia) — V ab i W b (x — ia))/2. The linear in A term now contains Ai(x) + Ai(x — ia), and 
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just gives the A field part of the continuum Df b = 5 ab di — gjah c A\. The quadratic in 
A terms perform the renormalization of <9j and give exactly twice the corresponding 
contribution to /tift 2 , because in that case only half of the expression arose from W 
fields which are parallel transported. Hence we find k 3 = Z^. 
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